{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we see how to call Gobnilp from with an R session and how to use Gobnilp and bnearn together. First of all let's attach the bnlearn package (which I will assume is already installed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: ‘bnlearn’\n",
      "\n",
      "The following object is masked from ‘package:stats’:\n",
      "\n",
      "    sigma\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(bnlearn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To communicate with Gobnilp we will use the excellent `reticulate` package (which is also assumed installed):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(reticulate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that reticulate is in use we can use the \"import\" command to import the `gobnilp` module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "gob <- import(\"gobnilp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the newly created `gob` to access attributes of the `gobnilp` module. Usually the class `Gobnilp`\n",
    "is all we are interested in and we just want that to create a `Gobnilp` object. Here's how to do that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "m <- gob$Gobnilp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that instead of Python's '.' to get to attributes we have R's '$'. Although in the above we created the variable `gob` usually there is no need to do so, since we just want the `Gobnilp` object, so doing this is simpler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "m <- import(\"gobnilp\")$Gobnilp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the object `m` we can use it to, for example, learn a BN from a data file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "m$use_discrete_data('discrete.dat',plot=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I have put `plot=FALSE` since although when running Gobnilp from R the learned BN is indeed plotted, errors are also generated. Note also that not even any textual output about the learned BN is here. This is an issue with Jupyter since this is textual ouptut *is* generated when using a normal R session. But we can still get a report on what was learned as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "**********\n",
       "BN has score -24028.09477835351\n",
       "**********\n",
       "A<- -5502.137377150637\n",
       "B<-A -3688.9395212202216\n",
       "C<- -3501.5105385969146\n",
       "D<-A,C -3555.0144442365527\n",
       "E<-B,F -4310.3049564706525\n",
       "F<- -3470.18794067853\n",
       "**********\n",
       "bnlearn modelstring = \n",
       "[A][B|A][C][D|C:A][E|F:B][F]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m$learned_bn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "**********\n",
       "CPDAG:\n",
       "Vertices: A,B,C,D,E,F\n",
       "A-B\n",
       "A->D\n",
       "B->E\n",
       "C->D\n",
       "F->E\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m$learned_cpdag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use Gobnilp to learn a BN from data supplied as an R data frame. Since we have bnlearn available we can grab its built in dataset `learning.test` and then have a look at it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " A        B        C        D        E        F       \n",
       " a:1668   a:2362   a:3717   a:1755   a:1941   a:2509  \n",
       " b:1670   b: 568   b:1024   b:1570   b:1493   b:2491  \n",
       " c:1662   c:2070   c: 259   c:1675   c:1566           "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data(learning.test)\n",
    "summary(learning.test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So `learning.test` is a data frame with 6 discrete variables. Here's how to use Gobnilp to learn from this data frame (with default parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2 <- gob$Gobnilp()\n",
    "m2$use_discrete_data(learning.test,plot=FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "**********\n",
       "BN has score -24028.094778353498\n",
       "**********\n",
       "A<- -5502.137377150637\n",
       "B<-A -3688.9395212202144\n",
       "C<- -3501.5105385969146\n",
       "D<-A,C -3555.014444236549\n",
       "E<-B,F -4310.3049564706525\n",
       "F<- -3470.18794067853\n",
       "**********\n",
       "bnlearn modelstring = \n",
       "[A][B|A][C][D|C:A][E|F:B][F]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m2$learned_bn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we get the same BN as before, which is no surpise since `learning.test` and the data in the file `discrete.dat` are exactly the same! Since we have only 6 variables there is no need to use Gobnilp's default parent set limit of 3. So let's learn again without it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2 <- gob$Gobnilp()\n",
    "m2$use_discrete_data(learning.test,plot=FALSE,palim=99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "**********\n",
       "BN has score -24028.094778353498\n",
       "**********\n",
       "A<- -5502.137377150637\n",
       "B<-A -3688.9395212202144\n",
       "C<- -3501.5105385969146\n",
       "D<-A,C -3555.014444236549\n",
       "E<-B,F -4310.3049564706525\n",
       "F<- -3470.18794067853\n",
       "**********\n",
       "bnlearn modelstring = \n",
       "[A][B|A][C][D|C:A][E|F:B][F]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m2$learned_bn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the same BN yet again. So this really is a globally optimal BN. We can use one of bnlearn's structure learning algorithms on `learning.test`. Here's hill-climbing, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "m3 <- hc(learning.test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "  Bayesian network learned via Score-based methods\n",
       "\n",
       "  model:\n",
       "   [A][C][F][B|A][D|A:C][E|B:F] \n",
       "  nodes:                                 6 \n",
       "  arcs:                                  5 \n",
       "    undirected arcs:                     0 \n",
       "    directed arcs:                       5 \n",
       "  average markov blanket size:           2.33 \n",
       "  average neighbourhood size:            1.67 \n",
       "  average branching factor:              0.83 \n",
       "\n",
       "  learning algorithm:                    Hill-Climbing \n",
       "  score:                                 BIC (disc.) \n",
       "  penalization coefficient:              4.258597 \n",
       "  tests used in the learning procedure:  40 \n",
       "  optimized:                             TRUE \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So hill-climbing with the BIC score finds the same network."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use Gobnilp to learn from continuous data in the same way. Here is an example using the bnlearn dataset `gaussian.test`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "data(gaussian.test)\n",
    "m4 <- gob$Gobnilp()\n",
    "m4$use_continuous_data(gaussian.test,plot=FALSE,palim=99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "**********\n",
       "BN has score -53258.94161814129\n",
       "**********\n",
       "A<- -7124.782936593145\n",
       "B<-D 487.26955605315015\n",
       "D<- -14692.560410628308\n",
       "C<-A,B -3743.043565645814\n",
       "E<- -10545.851006239516\n",
       "F<-A,D,E,G -7109.807736959381\n",
       "G<- -10530.165518128275\n",
       "**********\n",
       "bnlearn modelstring = \n",
       "[A][B|D][D][C|A:B][E][F|D:A:E:G][G]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m4$learned_bn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And again, we can see how hill-climbing does (using BGe scoring) ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "  Bayesian network learned via Score-based methods\n",
       "\n",
       "  model:\n",
       "   [A][B][E][G][C|A:B][D|B][F|A:D:E:G] \n",
       "  nodes:                                 7 \n",
       "  arcs:                                  7 \n",
       "    undirected arcs:                     0 \n",
       "    directed arcs:                       7 \n",
       "  average markov blanket size:           4.00 \n",
       "  average neighbourhood size:            2.00 \n",
       "  average branching factor:              1.00 \n",
       "\n",
       "  learning algorithm:                    Hill-Climbing \n",
       "  score:                                 Bayesian Gaussian (BGe) \n",
       "  graph prior:                           Uniform \n",
       "  imaginary sample size (normal):        1 \n",
       "  imaginary sample size (Wishart):       9 \n",
       "  tests used in the learning procedure:  75 \n",
       "  optimized:                             TRUE \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-53258.9416181405"
      ],
      "text/latex": [
       "-53258.9416181405"
      ],
      "text/markdown": [
       "-53258.9416181405"
      ],
      "text/plain": [
       "[1] -53258.94"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m5 <- hc(gaussian.test,score=\"bge\")\n",
    "m5\n",
    "score(m5,gaussian.test,type=\"bge\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hill-climbing and Gobnilp's exact algorithm find BNs in the same Makov equivalence class, and since BGe is *score-equivalent* they have the same score. In this particular example, the only advantage of exact learning is that we *know* we have the optimal BN."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
